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Abstract. A statistical inference method is developed and tested for pairwise 
interacting systems whose degrees of freedom are continuous angular variables, such as 
planar spins in magnetic systems or wave phases in optics and acoustics. We investigate 
systems with both deterministic and quenched disordered couplings on two extreme 
topologies: complete and sparse graphs. To match further applications in optics also 
complex couplings and external fields are considered and general inference formulas 
are derived for real and imaginary parts of Hermitian coupling matrices from real and 
imaginary parts of complex correlation functions. The whole procedure is, eventually, 
tested on numerically generated correlation functions and local magnetizations by 
means of Monte Carlo simulations. 
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Classical XY spin models with linearly interacting spins have been subject of 
intensive study in statistical mechanics since the mid 60’s concerning the investigation 
of critical phenomena [Tj. In particular, starting as classical lattice field proxy for 
the quantum theory of the A-transition of the Bose condensation of optical phonons 
[2], and further of liquid helium to its super fluid state [3], these models have been 
employed in the theoretical description of the 2D Kosterlitz-Thouless transition to an 
ordered unmagnetized spin vortex phase IS 15115] and, more generally, to the study 
of the chirality transition [3 El El [10] and to the critical behavior on random graphs, 
whose nature depends on their spectral dimension PM US- Further applications can 
be found to the roughening transition of the interface of a crystal in equilibrium with 
its vapor [[13] and to synchronization problems approached by means of the Kuramoto 
model [13 IH E], just to make a few examples. 

Pairwise XY models can, as well, describe propagation and amplification of linear 
waves in open and random media, as derived, e. g., in Refs. mmm as a classical 
degradation of the quantum theory for modes with overlapping resonances. In this 
framework, coupling constants yield information about the interaction between localized 
(inner) modes with discrete frequencies and radiative (outer) modes whose frequencies 
take values in a continuous dominion. In presence of relevant amount of disorder, the 
refractive index strongly and inhomogeneously depends on the spatial coordinates of 
the randomly placed scatterers inside the optically active medium. In these cases, a 
contribution to the couplings comes from the spatial overlap between the electromagnetic 
fields of the inner eigenmodes, modulated by a linear, inhomogeneous, susceptibility 
[211 [22]. A quantitative estimate of the coupling coefficients, thus, yield fundamental 
information about the space localization of the modes, and about the space dependence 
of the optical susceptibility. Eventually, XY pairwise models can be applied to another 
problem concerning light propagation, that is, the optimization of the transmission 
matrix of complex random media [23J [24], modeled as a coupling between input and 
output mode phases. 

In the present work we undergo the investigation of statistical inference techniques 
on XY models to provide a methodological theoretical frame straightforwardly applicable 
to the above mentioned problems with a particular focus on optics. The developed tools 
can be applied, as well, to any pairwise interacting model whose variables can be found 
in p states that can be considered as discrete values of an angle, the so-called p-clock 
model [251125]. Recent analysis has, indeed, shown that in the large (but not so large) 
p limit the XY model properties are promptly recovered for finite temperature and, 
further, very interesting features arise at finite small p [271 EH] . For p = 2, eventually, 
one recovers the Ising (i.e., Boolean) model for which inference studies have been carried 
out in Refs. [29l [30] . 

The paper is organized as follows: in Sec. [0 we first consider the inverse graphical 
problem [3lj for a real valued XY model. In Sec. [2j we present the study of the generic 
complex amplitude model with complex valued, and possibly disordered, interaction 
couplings and we derive the relationship between measurable correlation functions and 
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theoretical mode couplings under the hypothesis of a complete graph, i.e., in the so- 
called fully connected mean-held limit |29] were each spin is (feebly) coupled to all 
the others. In Sec. [3] we, further, test the obtained inference formulas on correlations 
numerically generated by means of Monte Carlo simulations on different kinds of 
underlying interacting networks, such as complete and sparse random graphs. We 
present our results on the efficiency of the proposed method both in the case of ordered 
exchange interaction and disordered couplings. Eventually, in Sec. 01 we discuss the case 
where the data sets for measuring the correlation functions, from which the couplings can 
be inferred, are small and in Sec. [5] we draw our conclusions and outline the perspectives 
of our work. 

1. Inference in XY model with real interaction couplings 

The simplest model we are going to consider consists of XY spins a = (cos (ft, sin (ft) with 
angles (ft 6 [0, 27t) , pairwise real valued interaction J l3 between sites i and j and an 
external held hi. Its Hamiltonian reads 

H = Jij COS (0j - (j>j) - ) h n COS (0i) (1) 

(t?) * 

where the set of interacting pair of sites (ij ) is determined by the topology of the 
underlying network and J rj is a symmetric matrix (the graph in undirected) whose 
elements can take any real value, deterministic or randomly distributed. The hrst 
focus of the present work is to derive the relationship between the two-point correlation 
function 

Cij = (<?i • &j) = (cos f>i cos (ftj + sin (fti sin (ftj) = (cos (cfti — cftf)) (2) 

and the interaction couplings Jij to infer the latter for the hrst. In Eq. (|2]) the average 
(...) is the ensemble average over the equilibrium distribution. Given the set of spin- 
spin correlation functions from experimental data, an inverse statistical problem is 
setup to investigate the interaction couplings among the spins in this model. Such 
inverse problems have been widely studied for finding: parameters of the discrete spin 
models using mean held theory for complete graphs |2!) [3(2 [32], 03], structural properties 
of proteins from multiple sequence alignment data, [33, 05], 36] effective local brain 
topologies from in-vivo neural recordings E7j. 

1.1. Variational free energy method approach 

Although the variational free energy method for XY model is somehow standard [138], 
we will briefly recall it here for fixing the notations. We aim at finding the probability 
distribution p((fti) by introducing a Lagrange multiplier A; to minimize the free energy 
subject to probability closure constraint JT / Q 27r d<ftip(<fti) = 1. The expressions for the 
internal energy E, entropy S and free energy F are the following: 

/•27T /»27T 

E = - / / d(ft i d(ft j p((fti)p((ftj)Ji j cos((fti-(ftj) 

i<j Jo J " 
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r-2n 


^ / d<f>ip((j)i)hi cos(0j) 


s=-E 


r-2 7T 


d^ip{(f>i) lnp(0j) , 


F = E-TS~y A, Qf 2 *, 


d(j>ip(<f>i ) 


(3) 


Taking the functional derivative of the free energy functional with respect to p(0j) and 
setting 5F/Sp(4>i) = 0, one finds the probability distribution p(4>i)- Let us first define: 

Hf = y Jjj (cos((j6j)) + ^ , II 1 / = y J i:j (sm((f)j)) , 


H % = J(Hff + (Hf) 


y \2 


a, = arctan ■ 


H[ 

HI 


and the normalisation factor: 

/■2ir 


Z = 


d0;exp (Hi cos (fi - a/)) = I 0 (Hi ) 


(4) 


(5) 


'0 


where J 0 (T) is the modified Bessel’s function of the first kind. Substituting Eqs. (j3J) and 
(J5} in 5F/5p((j)i ) = 0 one obtains 

exp ( HiCos((j)i - cti)) 


p{4>%) = 


Io(Hi 


( 6 ) 


The magnetization components are, then, derived averaging cos(0 7 ) and sin(^j) over 
the probability measure p(<j>i), yielding 

< = (cost*)) = , m? = <*.(*)) = IBM ( 7 ) 

-'ol-'uj H [Hi) 

where the modihed Bessel’s function of the first kind and their derivatives read: 

r 2n 


h(z) = / #cos(0)e 2COsW 

Jo 

Iq(z) = h(z) , l[(z) = I 0 (z ) 


( 8 ) 


AW 


Correlation functions are consequently computed using the linear response formulas 
i. e., deriving magnetizations with respect to perturbations in the external fields h k and 
yielding 


C X _ 

ik ~ 


8m? 


h(Hf)Hl 


8h k 8h k \ I 0 (Hi) H 

Hf fl\_I[\ SH t | h SHf Hf /, SH 
r o II 


+ 


Hi \I 0 /q / 5h k IqH, Sh k Hf 7 0 Sh k 
(H x )f , , s(H y )r 


S_H[ 

8h k 


0 q% - \m\)- 


Hf v ‘ Hf 


SHVHfHV , n _ 
8h k Hf 




+ (W) : 

2 -i) 


( 9 ) 
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r y _ 5 rrv[ 

ik Sh k 


6H y 


(ft 


m. 


(-H 


y '|2 


Sh k 

SHf hthj 


' ' Hf 

(2ft - 


Sh k H\ 

where we make use of the following substitutions: 

VI ! i w , h 


m. 


+ ti¬ 
ll 2 ~l) 


ft) 


(H 


x\2 


m 


qi = (cos 2 (0i)) = 1 


In Hi ’ 


mi \ 2 = (m x ) 2 + (m y ) 2 


We further dehne 


m 

hi = ^ 
ml 


I |2 7 

\ m i\ = ~r 

to 


y 


! 1 

n 


= cos a* = cos[arctan(/ij)] = 


HI 
H 
H y 

-j- = sin cq = sin[arctan(/ij)] = 
Hi 



fl\qii m i) = 
f2 l \qu m i ) = 


\/l + hi 
mf- + p 2 (l - q^ 


1 + hi 

(1 + /i 2 )(l - ft) + ^(1 - |m;| 2 ) 


g M (qi,mi) = v ^ 


2ft 


m. 


1 + h* 2 
2 -l 


1 + hi 


( 10 ) 


( 11 ) 

( 12 ) 

(13) 

(14) 

(15) 

(16) 

(17) 

(18) 


Plugging the derivatives of Hf and Hf with respect to h k into Eqs. fl9l)-(fTT)l) and making 
use of Eqs. dill) - (fT8j) . we obtain a system equation for the correlation functions C x f y 
and the couplings matrix J lk . To be compact we adopt the operator form for expressing 
observables equivalent to co-ordinate form: J is the couplings matrix, C x,y the two point 
correlation matrices, I is the N x N identity matrix and for any vector z = {z\, ..., zjv}, 
the N x N matrix 1? is the diag(z) matrix, i.e. a matrix with the elements of z on the 
diagonal and all other elements equal to zero. The system ©-(HOD, thus, reads 

C x = I fi [JC X + I] + lgJC y (19) 

C v = % [ JC X + I] + I f - 2 JC y (20) 

Let us consider the substitutions 


1a = JC x + l (21) 

1b = JC y (22) 

Eqs (fl9l) and (120T) for the correlation functions can be, accordingly, rewritten as 

C x = + lglB (23) 

C y = + Iffts (24) 
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To infer interaction couplings of the system we solve equations ff2Ul) - (EHl) in y’s yielding 


lA = h c*-i % c y 

1 B = \C X - i k c v 


where 



A, = fPf? ~ (9 (<) ) 2 


We, eventually, obtain the values of J by inverting Eq. (121)1 

J=(7a-I)(C")- 1 



(25) 

(26) 

(27) 

(28) 

(29) 


and using 7 ^ obtained from measured correlations and magnetizations. Substituting for 
Jij in Eq. ([4]) and using Eqs. (jT4lb(lT5l) we, moreover, obtain the inference formula for 
the external field from the inferred J’s and the measured m’s: 


Iph = Jm y - I p Jm x 


(30) 


1.1.1. Zero external field Considering the case at h 
simplify as 






H? = 0 


a,- = 0 


x _ h{Hi) 

1711 /o (HiY 

Ai ) 2 

/l — Qi ~ m i 


, m y = Hi = 0 

Qi, 


f? = 1 


9i = 0 


0 , Eqs. ©, dZD, (HID-(USD 
Hi = Hf (31) 


and the correlation function reduces, then, to 
Srnf 


C X _ 

ik ~ 


Sh k 


h=0 


[(cOS 2 (0i)) - (COS(^)) 2 ] ^ J il C jk + 5 ik 


Further, inverting Eq. (1321) . eventually yields: 


(32) 


(33) 


In the above Eq. (1331) of inferred J, it is worth noticing the similarity with the expression 
of inferred J in terms of C of the Ising model with discrete spins, see, e.g., 
where (cos 2 (f>i)) = 1 . 


2. Model with complex spins and couplings 

In this section, the interaction coupling matrix J tj is considered to be a complex matrix 
consisting of a real and an imaginary part as = J[j + i./f and the external held is 
a complex vector hi = hf + ih\. This model can be a proxy for the propagation and 
interaction of waves in optically active media, ordered or random, where electromagnetic 
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modes can be represented, in a properly defined base of eigenvectors, by complex 
numbers a* = e 1 ^: each mode amplitude has its own magnitude A , = \ai\ and 

phase angle 0* = arg(aj). Indeed, the electromagnetic field can be decomposed, in the 
slow amplitude approximation [39], in terms of the complex amplitudes of the modes 
localized inside the medium 

N 

E(r , t) = J2 a k (t) U k (r) e l + c.c. (34) 

k=1 

where the frequencies oj k take values on a discrete dominion and U k s are the eigenvectors 
of the eigenmodes in some given basis allowing for a decomposition between inner and 
outer modes by means of Feshbach projectors [20]. We stress that, in cavities with non- 
negligible leakages or cavity-less light scattering random media, the N modes indicated 
in the sum are by no means a complete basis, but they are the subset made of purely 
localized modes amplified inside the cavity. 

Such modes can be proven to display a stochastic dynamics governed by a quantum 
Langevin dynamics [T9] , In the classical limit such evolution is proved equivalent to the 
master equation for the density of states ra and in terms of complex amplitudes takes 
the form 


dn(t ) ^ ' Uinvflrri: 4 “ Qq* 

m m 

{rjn(t )) = 0, (Vn(f) r ln'(t')) ~ 2 n t hS n md(t - t') 

where, for large enough heat-bath temperature, the thermal number of photons in the 
classical regime is n t h oc T and the complex valued white noise is approximated as 
uncorrelated on different states [40, H]. In general, though, we recall that non-diagonal 
covariances in the space mode can be non-zero for Markovian dynamics GDI- The linear 
non-diagonal coupling J nm between modes is the dumping matrix, associated to the 
openness of the optical cavity, due, e. g., to leakages in standard lasers 02! or to the 
cavity-less structure of the scattering region of the optically active material in random 
lasers [43( [44] [45] [21] . The total power taken by the system is a constant £ that, rescaling 
the amplitudes as a n a n /y/cJf, can be expressed as a simple spherical constraint 

N 

£ = eN = |a n | 2 (36) 

n =1 

The static properties of the above dynamics can be derived by studying the 
Hamiltonian 

^[W] = -Y. a fhjO* - ^2 h,a* (37) 

(y) * 

where Jij includes the dumping incorporating the inner-outer modes interaction, 
and, possibly, also includes the spatial overlap of the eigenmodes modulated by an 
inhomogeneous dielectric constant. We add a complex external field hi for generality. 
In the present work, we assume J tj to be Hermitian, thus = J and Jh = — jl. 
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The dynamics of mode phases changes at a much faster time scale than mode 
amplitudes [36] . Moreover, in presence of a large number of modes and a not too diluted 
interaction network, in a wide variety of systems intensity equipartition occurs among 
all modes, i.e., \a t \ ~ 1 m, trivially satisfying constraint (136|) . In particular, pairwise 
interacting mode systems, in any graph topology, always display intensity equipartition 
in all thermodynamic phases. One can thus work in the so-called quenched amplitude 
approximation, that is, amplitudes are further taken as quenched and incorporated in 
the J’s, yielding the Hamiltonian: 


H[{())}\ = - Y i J S cos (^ “ fa) + 4 sin (0* “ fa)] 

(ij) 

~ Y \- h i cos (^) + h i sin(0i)] (38) 

i 

Though derived in terms of light modes interacting in an optical open cavity, we 
stress that the above Hamiltonian generically describes any linear wave system. Indeed, 
it can, e. g., describe a class of optimization problems where the variable 0 represents 
the phase of a pixel of the incoming/outcoming light propagating through a random 
medium, including disordered optical fibers, and -fj represents the transmission matrix 

mm- 

The same variational procedure of Sec. [His applied to the complex system modeled 
by Eq. f[38]h The following substitutions for Hf and Hf are considered, in place of Eq. 
(J3j), in order to calculate the probability distribution p(<f>i): 

Hf = Y J ij( cos (fa)) ~ Y 4< sin ( < 4> + h * ( 39 ) 

j i 

H f =Y J 5( sin (fa)) + Y 4< cos ( < 4> + h i ( 4 °) 

j j 

It is found that the structure of p(0j) remains the same as in the case of real-valued J’s, 
cf. Eq. (jnj), though the components Hf,Hf and Hi adorn different expressions. 


2.1. Correlation functions and inference formulas 


In this section, equations for the correlation functions are derived. In the presence 
of complex fields, four correlation functions are found from the differentiation of 
magnetizations with respect to both components of the fields externally acting on the 
system. Using the expressions of Hf /Hi and Hf /Hi as in eqs. (TT4T) - (TT5T) . magnetizations 
in Eq. (J7J) can be written in the following form, 


m, = 


W<)h; 


mf = 


h(Hi)H? 


Io(Hi) Hi ’ 8 I 0 (Hi) Hi 

To derive correlation functions we use the following linear response relations, 


(41) 


r x _ forcf 

ik Sh£ 


C x 

ik — 


dmf 


v _ 5m%_ 
ik 6h «’ 


Ay _ Yl 
ik 5 hi 


( 42 ) 
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Performing the above derivatives as in Sec. [U we find equations for the correlation 
functions in matrix form as following 

C x = I fi [J R C X - J ! C y + I] + % [j R C y + J : C X ] (43) 

C x = I f - \j R C x - JcA + % \j R C y + J I C X + l] (44) 

C y = % [J R C X - J ! C y + I] + I f - [J R C V + J J C X ] (45) 

C y = % \j R C x - J J C V 1 + \j R C y + J J C X + l] (46) 

These are decoupled two-by-two and the two subsystems are not independent, but 


equivalent to each other. Therefore, to obtain the inference formulas for J we can 
simply solve Eqs. (jj3j) and (l45]h To invert them let us first define: 


r 4 = J R C X - J 7 C y + I (47) 

r B = J R C y + J ! C X (48) 

Substituting into eqs. (l43lh (1451) and solving for T’s we obtain 

r.4 = \C X - \C y (49) 

T B = \C X - \c v (50) 


where coefficients are given in Eq. (l27]h Now, after obtaining T^ and T^ in terms 
of measurable quantities, we get back to eqs. (1471) and (|48li and solve them to extract 
interaction couplings, yielding the main equations of our work: 

j R = [(r A -1) (c v y l + v B (c x y l ] [c x ( c y y l + c y {c x y 1 )} ~ l (51) 
j 1 = [- ( r A + 1 ) (c 1 )" 1 + {cy- 1 } [c x (c y ) _1 + c y (c x ) -1 ] _1 (52) 

Eventually, to infer the external field values from the inferred J values and the 
measured magnetizations, using H x = H y /^ in Eqs. (|39li - ()4(l . we obtain 

I- h R - h 1 = ( 1 m x + % myj 1 (53) 

3. Numerical tests 

To verify the efficiency of the predictions of the inference method derived above, we 
present our tests on data provided by means of Monte Carlo simulations of models 
exactly given by Eqs. (pQ) and (1381) with known couplings and fields and defined on 
different topologies. 

We will first consider data generated from finite size models on fully connected 
long-range models, where each spin is connected to all the others and each coupling J 
is vanishingly small to guarantee thermodynamic convergence. These are the hypothesis 
at the ground of the variational approach followed in deriving Eqs. fl51H52p . We will 
consider both the case of ordered Jij = J and randomly distributed Jij with Gaussian 
probability of mean zero and variance one. 
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We will afterwards consider data generated on models defined on Erdos-Renyi 
(ER) sparse graphs [48] where the connectivity is randomly distributed according to 
the Poisson Distribution 

m = — ( 54 ) 

Once again, both deterministic and quenched disordered interaction couplings are 
considered. 

The values of the inferred matrix are, eventually, compared to those of the original 
matrix. All the data, namely, correlation functions C X ,C Y , magnetizations m x ,m y 
and q and \m \ 2 used for this analysis are computed from therm a h zed data produced 
using Monte Carlo simulations with the parallel tempering algorithm. For each case, we 
display (i) the comparison of the whole range of sorted original and inferred couplings, 
(ii) the comparison of inferred to original couplings to one single site and (iii) the 
sensitivity plots for true positive (fraction of original non-zero couplings inferred to 
be non-zero) against the number of predicted connections, yielding an insight into the 
topology of the graph. 

3.1. Ordered couplings on complete graph 

For the fully connected case the entire analysis is shown at T = 0.5, in the ferromagnetic 
phase. In Fig. |T|we display the comparison between the original (dotted/green lines) 
and the inferred (continuous/red and black lines) couplings sorted by magnitude in a 
system of N = 64 modes and N(N — l)/2 = 2016 independent couplings. We consider 
three cases. On the left panel purely real, ferromagnetic Jij = Jji = J couplings are 
plotted in a zero external magnetic held, h — 0, and for h = 0.2. To infer the values of 
Jij we used Eq. (l5Tf that fairly predicts the initial couplings even though the equations 
used are of the general case, i.e., complex Hermitian J’s that can take any value. Using, 
instead, Eq. (1291) or Eq. (j33j) . focused on the specific cases of purely real couplings and 
zero-held, no difference is appreciated down to the third digit. 

In the mid and right panels, we plot real and imaginary inferred J’s with zero helds 
and for h R = 0.4, hi = 0.2. The quality of the prediction is comparable with the purely 



Figure 1: Sensitivity plot of the sorted inferred interaction couplings in the complete 
graph at T = 0.5 for an ordered system of N = 64 variables. Left: Ordered, only real J 
for h — 0, 0.2. Right: J R and J/ for h = (h R , hi) = (0, 0) and for h = (0.4, 0.2). 
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Figure 2: Sensitivity plot of the sorted inferred interaction couplings in the complete 
graph with IV = 64 and quenched disorder at T = 0.5, in the spin-glass phase. Left: 
Purely real valued J for h = 0,0.2. Right: Real and imaginary parts of the complex 
disordered couplings, Jr and Ji for h = (hR,hj) = (0.4, 0.2). Bottom panels: original 
and inferred couplings to a single site. 


real case. 

3.2. Disordered couplings on complete graph 

We, then, inferred couplings from correlation functions and magnetizations generated 
in a system of N = 64 spins where Jr and Jf s are originally generated by means of a 
Gaussian random distribution of mean zero and variance one. In the left panels (top 
and bottom) of Fig. [2] the case of purely real couplings is exposed, both in h = 0 and 
h = 0.2. No difference is appreciated between these two cases and the magnitudes of 
both are about the magnitude of the original couplings. As detailed in the bottom panel 
for couplings to a specific site, inferred J’s faithfully predict sign and magnitude of the 
original ones. The center and right panels display the behavior of real and imaginary 
part of a system with complex couplings both in absence and presence of external fields. 
Again, the presence of external fields do not alter the inference predictions and signs 
and magnitude of original couplings are correctly predicted. 

3.3. Ordered couplings on sparse random graph 

Next, we show the analysis for the case where the J matrix is diluted, though the 
formalism developed in this work is rigorous for fully connected systems and not for 
sparsely connected ER graphs. The connectivity probability is randomly distributed 
according to Eq. (I54|) with average connectivity c = 6. Data shown are for N = 256 
spins at temperature T = 2.5 > T c ~ 1.95 for the systems with complex Hermitian 
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Figure 3: Sensitivity plot for the inferred couplings of an original random ER graph of 
N = 256 nodes with fixed deterministic values of J at T = 0.75. Left: Real valued J’s 
for the largest entries (first N c/2 elements, c = 6). Right: Real and imaginary part of 
complex couplings (first 3/2N c entries). Bottom panels: original and inferred couplings 
to a single site. 


couplings, cf. Eq. (l38]l and T = 0.75 (here T c = 2.958) for systems with purely real 
couplings, given by Eq. ([I]). 

In Fig. [3] we display the comparison of the inferred and the original J’s by means of 
Eq. (f5TT) . In the left panels original couplings are all real and, when non-zero, all equal 
to each other. Each site is connected to a finite, At-independent, number of others, c = 6 
in the average. The analysis gives correct indication for non-zero Jr s both in absence 
and presence of external fields. However, our method always provide non-zero (though 
small) predictions for all couplings. Indeed, the true positive plots, cf. Fig. [4| decay 
down to zero only gradually, not sharply, quantifying the wrong predictions.Even though 
all non-zero elements of the matrix have been predicted correctly, for every zero element 
the formalism does not predict exact zero, bringing down the score of true positive. 

The same situation arises for complex couplings, where rather good estimates of 
non-zero J R and J 1 entries is provided, including a sharp decrease of the value of the 
inferred couplings at the sorted coupling N c/2, cf. top panels in Fig. [3] This is, 
though, contrasted by the rather poor estimate of zero couplings. As confirmed by the 
true positive plots in Fig. [4l the right panels of Fig. [3] show that zero couplings, i. 
e., those beyond the N c/2-th coupling, are inferred to acquire a non-zero value. Zero 
original couplings are not reproduced at all in the sparse case. 
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Figure 4: True positive plots for the ER sparse random graph with average connectivity 
c = 6, size N = 256 and with deterministic couplings J — 1. 


3-4- Disordered couplings on sparse random graph 

The inference maintains the same quality also in the case of random values of the 
couplings. In Fig. [5l left panels, we display the case of a ER random graph whose 
couplings have Gaussian distributed real values, with average zero and variance equal 
to one. The top figure in the sensitivity plot for the first N c/2 couplings, with 
and without external field, compared to the original disordered coupling values. The 
bottom panel show the comparison between original and inferred couplings to the graph 
node 3, to exemplify that: (i) all original non-zero couplings are well reproduced and 
discriminated in the inference procedure and (ii) all inferred couplings are non zero, also 
those corresponding to missing original couplings, though the latter acquire a rather 
small value in comparison to the inferred true links. The same analysis is illustrated 
in the mid and right panels of Fig. [5] for the real and imaginary part of a system with 
Hermitian couplings. In the sensitivity plot the first and last 3 N c/2 couplings are 
reported and compared to the original ones, signaling that the inference quality is very 
good, though non-zero couplings are inferred to have a small non-zero value. In the 
bottom panels couplings to one node are displayed. This is confirmed in Fig. [6] where 
the true positive curve is shown to decrease sharply after the last non-zero coupling but 
still is non-zero for all matrix entries in all considered cases. 
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Figure 5: Sensitivity plot for the inferred couplings of an original random ER graph 
with Gaussian distributed random couplings (zero mean, unitary variance) at T = 2.5 
and with N = 256. Left: Real valued J’s for the largest entries (first N c/2 elements). 
Right: Real and imaginary part of complex couplings (first and last 3/41V c entries). 
Bottom panels: original and inferred couplings to a single site. 





Figure 6: True positive plots for the ER sparse random graph with average connectivity 
c = 6, size N = 256 and with Gaussian random distributed couplings (zero mean, unit 
variance). Left: True positive of a system with purely real couplings and zero external 
held. Mid and Right: True positive of the real and imaginary parts of the couplings in 
a system with zero held. 


4. Small data size behavior 

In this part, we show how the quality of inference is deteriorated as the number of 
measurements composing the data set used to calculate correlations decreases. In 
the main hgure 0 the entire sorted J matrix is shown and in the inset the absolute 
value of the hrst 2000 elements are shown. We see that the sensitivity plot remains 
the qualitatively the same for all data sets, but the transition from non-zero to zero 
couplings becomes sharper and sharper as the data size increases, yielding evidence for 
an underlying sparse graph. 
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Figure 7: Sensitivity plot of data generated from an ER graph of N = 256 sites and 
mean connectivity c = 6 with real ordered couplings and no field. The behavior of the 
original network is marked as a full black step curve. Correlation functions averaged 
over, respectively, 4096, 8192, 32768 and 131072 data from Monte Carlo simulations are 
considered. In the inset a detail of the absolute values of the bonds around the last 
non-zero original coupling (N c/2 = 768) is provided. 


5. Conclusions 

In the present paper we have derived an inference procedure to determine the 
coupling constants of pairwise interacting systems with continuous XY spins, complex 
interactions and complex external fields. 

For testing the analytic inference approach we have applied it to data numerically 
generated by means of Monte Carlo simulations at equilibrium and we have compared 
the inferred values of the coupling constants to the ones of the simulated system. We 
considered models with disorder in the coupling values and models with disorder in 
the coupling connectivity, studying both complete and sparse random graphs with both 
deterministic and quenched disordered couplings. The inferred couplings turn out to 
reproduce original ones in an excellent way in fully connected models, that is, under the 
conditions at the ground of the theoretical derivation of inference formulas Eqs. ()29l I5T1 
l52]h Also when applied to sparse random graphs, though, the quality of the inference is 
of a high standard. The only problem arises in the wrong representation for missing links 
that acquire always non-zero value in the inference procedure. Their values, actually, 
decrease with increasing data size, but do not reach zero even for very large data sizes, cf. 
Fig. [7] Else said, the true positive curve is always non-zero for all couplings. The value 
of false positive inferred couplings turns out, though, to be systematically much smaller 
than the value of true positive bonds, with a sharp quantitative distinction between the 
two. 
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In the field of random photonics, the reported method to quantitatively infer 
coupling constants from experimental data would allow to obtain estimates of the 
effective damping interaction between localized modes in a random medium mediated 
by radiative modes [18] [Hi 20] in an open cavity, to extrapolate the magnitude of the 
optical-response-modulated spatial overlap between those modes and, eventually to 
obtain information about their localizations. Further on, the inference method for XY 
pairwise models can be applied to the optimization of the output signal from complex 
random media [23., [24], including disordered optical fibers, by inferring the elements of 
the transmission matrix. 

Further investigation on inference of waves can include alternative probes of the 
linear problem here reported by means of other inference methods such as pseudo- 
likelyhood. Most interesting is the generalization to nonlinear problems, allowing for 
the reconstruction of the properties of light modes in both ordered and random lasers, 
both in the continuous and in the pulsed regime [27] [47]. 
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